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Brief Description 

A FORTRAN CODE has been developed which predicts pressures and flow rates in a 
piping network for the regimes of laminar, transition, and turbulent flow. A novel and 
effective scheme for modeling the transition region between laminar and turbulent flow 
is demonstrated. A proposed linear indicator for the deviation from laminar behavior 
is introduced into the analysis in order to facilitate the computational task. Although 
the turbulent flow problem is inherently nonlinear, the method of analysis described 
herein, typically employed with linear electrical networks, is found to converge rapidly 
to an approximate solution without user- specified a priori guesses. 
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DUDLEY KNOX LJeRARY 
NAVAL POSTGRADUATE SCHOOL 
MONTEREY CA 93943-5101 

1. Introduction 

An interest in the properties of incompressible fluid flow in pipes can undoubtedly be 
traced as far back as the era of the Roman empire. A brief synopsis of the historical 
development of fluid mechanics can be found in many of the introductory level texts 
[Vennard (1961) for example]. 

The problem of accurately solving for the parameters in incompressible fluid flow 
in piping networks is dependent upon the single-pipe model for fluid flow. Many 
semi-empirical methods for characterizing the flow properties have appeared in the 
literature [Matthew (1981), Wuori (1993), and Churchill (1977)]. Since the early work 
by Feigenbaum (1978), see for example [Hofstadter (1981)], on chaos in simple sys- 
tems, considerable progress has been made on analytically characterizing turbulence 
in fluids [Landau and Lifshitz (1987)]. Nevertheless, many engineering problems are 
solved using the Stanton diagram, which shows Nikuradse’s measured data. 

For engineering problems which require repetitive numerical evaluation, an empir- 
ical relationship facihtates the task and various investigators [Tsang and Kee (1987), 
Round (1980, 1985), and Colebrook (1939)] have proposed such formulae. These for- 
mulae have the drawback that the transition region between laminar and turbulent 
flow is inadequately represented. Only the fluid flow properties in the laminar region 
are fairly well understood [Vennard (1961)] using methods of fluid analysis. 

One of the earliest and most well known methods for solving pipeline networks 
is attributed to Hardy Cross [Giles (1962) and Lindeburg (1982)]. The Hardy Cross 
scheme, not only requires an a priori guess on the initial flows, but upon comparison 
with more recently reported methods [Carnahan and Christensen (1972), Bending 
and Hutchinson (1973), Gay and Preece (1975, 1977)], it is found to be less efficient. 
Finally, in a recent Naval Postgraduate School (NPS) thesis [Ellis (1993)], a pipeline 
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method formalism based on the Cholesky method for matrix reduction heis been 
investigated. Although the mathematical analysis of the scheme is fairly efficient, the 
physical modeling did not encompass conditions for other than laminar flows. 

The main purpose of this report is to extend the range of flow to include transition 
and turbulent flow regimes. As this is primarily a feasibility study, various facilitating 
assumptions were made. First, it is assumed that English units are preferred for 
input/output. Second, all pipes are assumed to have uniform circular cross-sections. 

2. Physical Modeling 

For incompressible fluids, two of Bernoulli’s rules [Lindeburg (1982)] for conservation 
of mass and energy are, the continuity equation 

Q = A,V^=A2V2 (1) 



and Bernoulli’s law 



^ + (P2 - Pi) + pg{z2 - -2i) + pgh' 



( 2 ) 



Q " 2 

where the subscript 2/1 refers to positions 2/1 and the symbol lists given in Section 
10 defines the other relevant parameters. The head loss, as defined by the Darcy- 
Weisbach equation, is given by 



h' = 



fLV^ 

2gd 



( 3 ) 



where L is the pipe length, d is the hydraulic diameter, and / is the friction factor. 
For circular cross-section pipes the hydraulic diameter and the physical diameter are 
identical. Using similitude analysis [Vennard (1961)], it is possible to show that the 
friction factor depends only on two dimensionless parameters, the Reynolds number 

pVd 



Re = 



( 4 ) 



2 



where fi is the dynamic viscosity, and the relative roughness 



c = 



e 

d 



( 5 ) 



where e is the roughness of the pipe. 

A sourceless section of pipe which is horizontal {z 2 = 2 : 1 ), and with constant 
cross-section so that V 2 = Vj, will, from eq (2) and eq (3), then have a pressure drop 



or after use of eq (1) 



Pi - P 2 = h = pgh' = 






pfLV^ 

2d 



( 6 ) 



(7) 



This suggests the useful electrical analogy with Ohm’s law where pressure drop cor- 
responds to voltage and fluid flow corresponds to current. It follows that consistent 
with the analogy, the effective resistance of the pipe is then given by 

.LpV 



R = f 



2dA 



( 8 ) 



where the units of this resistance are defined by the ratio of pressure to flow. 

It is observed that the problem of predicting fluid flow in pipes now reduces to 
accurate characterization of the friction factor. For the laminar flow domain the 
analytically derivable [Vennard (1961)] result 






Re 



(9) 



agrees very well with measurement out to Re < 2100 which is the well known limit in 
pipes, ducts, tubes, and channels for laminar flow. After substitution of eq (9) into 
eq (8) it is found that: 



Rt 



32Lp 

<PA 



( 10 ) 
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and the corresponding relationship for head loss is 



ht 



d2Lfi 
d^A ^ 



( 11 ) 



This is in agreement with an expression cited by Bending and Hutchison (1993). The 
subscript £ has been used to indicate the laminar flow domain. Note from eq (10) 
that, as expected, the pipe resistance is independent of the flow rate in the laminar 
regime. 

In order to easily characterize the degree to which the fluid flow is turbulent, the 
head loss can be expressed as 

h = h( X xp (12) 



where the turbulence factor, denoted xp,- is a Hnear gauge of the deviation of the 
system from laminar flow. After noting that eq (7) can be reexpressed as 



H = X 



(13) 



d'^A '' 64/i 

it can be concluded from eq (4) that, in agreement with the definition of the Reynolds 
number, the turbulence factors satisfy 



<p = f- 

^ ^ 64 



(14) 



As seen from eq (9) ^ = 1 for laminar flow. In the next section it will be shown that 
for transition and turbulent flow ^ > 1. 

The pipe resistance cam now be expressed as 



R = Rt X xp 



( 15 ) 



where as previously noted, Rt does not depend on the flow rate. 
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3. Friction Factor Estimation for Nonlaminar 
Fluid Flow 

For purposes of the development here it is necessary to introduce two distinct non- 
laminar regions. First, the transition region satisfying 



Reo > Re > 2100 



(16) 



where Rep is usually cited typically between 3500 and 4000 [Vennard (1961)] and for 
the turbulent region 

Re>Reo (17) 



One of the most commonly employed empirical formulas for the turbulent flow domain 
defined by eq (17) is attributed to Colebrook (1939) 

1 



V7 



= -21og 



e 2.51 
3Jd Rev/J 



(18) 



Although, as seen in eq (18), the dependence of the friction factor on Reynolds number 
is implicit, the expression will converge rapidly upon iterative substitution. 

To date the most viable predictors for the trainsition region, see eq (16), are 
empirical. The Stanton diagram [Vennard (1961)] shows smoothly varying curves 
which are based on measurements in the transition region. The Moody diagram, 
which only plots expressions for eq (9), for Re < 2100, and eq (18), for Re > Re,,, is 
not useful for the transition region. Bending and Hutchinson (1973) suggested, in a 
computer analysis of piping networks, that a line 2 ir interpolation between the laminar 
and the turbulent domain be made to provide friction data in the transition region. 
For the network analysis scheme to be described herein, a linear approximation for 
the transition region was found to be insufficiently accurate upon numerical testing. 
In particular, the iterative program with the finear approximation would either not 
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converge or converge very slowly. The problem in this scheme was the dramatic 
discontinuity in the first derivative at the edges of the transition region. A greatly 
improved modeling scheme, in terms of convergence, is given by 



Re 

’ ~ 2T00 



[/(Reo) 



Re 



2100 J 



sin 



7T (Re-2100) 



2 (Re„ -2100) 



(19) 



where /(Re^) was calculated via the Colebrook model eq (18). Note that, hke the 
linear approximation model, continuity is preserved at the edges of the transition 
region. Because, as seen from eq (19), the slope at the transition edges is zero, the 
discontinuity in the slope is much less dramatic and the convergence was found to be 
much more rapid. Figure 1 is a graphical representation for these curves in standard 
log-log format for various relative roughness factors. 

It is worth noting that the work by Churchill (1977) agreed very well with the 
scheme just described for Re > 2100 and provided an explicit formula for the friction 
factor, /. However, for Re < 2100 (laminar flow) the expression deviated significantly 
from the classical expression of eq (9). Hence it was not used. 

The argument supporting the claim that the turbulence factor eq (14) is generally 
greater than unity is loosely based on the following observations. First, in 1913 Blasius 
proposed an empirical relationship for smooth pipe, i.e., c = 0 in the turbulent flow 
regime [Vennard (1961)] 

0.316 



yinnooth pipe — 






3000 < Re < 100,000 



( 20 ) 



Second, the construction algorithm for predicting / in the transition region, eq (19) 
requires that 



/itr&nsition ^ /uminar — 



( 21 ) 



This leads to the inequality that is estabflshed from a comparison of eq (20), eq (21), 
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and Fig 1 



64 

/rough pipe — /smooth pipe ^ /laminar 



( 22 ) 



and is supported by analysis and measurement [Vennard (1961)]. The inequality, 
^ > 1, then follows directly from the definition of the turbulence factor, given by 
eq (14). 

The turbulence factor hais been introduced into this report for two reasons. Al- 
though the friction factor is a linear indicator for the system response of the fluid 
dynamics, it is not a linear indicator for the deviation of the system from laminar 
behavior. Note that the quantity (Re — 2100) is a nonlinear indicator for the devia- 
tion from laminar behavior. The purpose of this work has been to extend the Ellis 
(1993) effort to encompass nonlaminar regimes. The turbulence factor will readily 
demonstrate the need for this extension. Second, it turns out that the laminar pipe 
resistance, i?<, can be computed at the onset of the problem and only the turbulence 
factors need to be updated in the iterative calculations of eq (15) rather than to com- 
pute eq (8). In brief, the turbulence factors provide some computational advantages 
which, for large piping networks, would be significant. 



4. Basic Circuit Modeling 

Following the electrical analogy for an arbitrary pipe, a pressure difference, Apk, 
across a pressure source, Ap^k, in series with a pipe resistance, rjt, with a flow qk, 
satisfies 

Apk = qkfk — Aptk (23) 

consistent with the convention established by Fig 2. The corresponding relation in 
terms of a current source is given by 



qk = q,k + VkApk 



(24) 
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Figure 1: Friction factor model. 1 
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Figure 2: Pressure source circuit model. 

where ?/* = 1/rjt. Figures 2 and 3 define the nomenclature and the topology. Equation 
(24) can be compactly represented for adl branches or pipes in matrix format 

Q = Q, +YAP (25) 

where it is seen that eq (24) and eq (25) have unresolved unknowns associated with 
flow rate and pressure. The condition of continuity at the nodes 

^9*, = 0 ; = l,2,...nT (26) 

* 

where nj is the total number of nodes in the network, provides enough information to 
resolve the flow rates and pressure drops. A formulation for this process is presented 
in the next section. 
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Figure 3: Flow source circuit model. 

5. Mathematical Formulation 

The mathematical formalism is presented via an illustration taken from the Ellis 
(1993) thesis. The piping network shown in Fig 4a has five nodes and six connecting 
branches. It is also seen that the branch between nodes 1 cind 5 contains a pump. 
This pump provides a source for fluid flow at its high pressure side shown at node-1 
with am arrow. The low pressure side is the corresponding sink. The corresponding 
flow graph is shown in Fig 4b. Each branch now has an arbitrarily chosen orientation 
arrow which defines the convention for positive fluid flow. Note that consistent with 
the flow direction shown for branch 1, the source shown on Fig 4a would be cissigned 
a negative numerical vade. The first step in the mathematicad formahsm is to cre- 
ate both a network representation for the fluid pipehne network and a corresponding 
flow graph. The flow rates in the b branches, in this case 6, can be represented as a 6x 1 
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(b) 



Figure 4: a) Six branch pipeline network; b) Associated flow graph for a) [15]. 



column vector 



Q = 



91 

92 

93 

94 



(27) 



98 



The elements of the augmented node-branch incidence matrix, c? , can be defined 



as 



= 


1 


Branch j leaves node i 


(28a) 


4 - = 


-1 


Branch j enters node i 


(28b) 


4 = 


0 


otherwise 


(28c) 



where i = 1 , 2, . . . , nj nodes and j = 1, 2, . . . , 6 branches. This leads to: 







1 


1 


0 


0 


0 


0 ■ 




• 




0 


-1 


1 


1 


0 


0 






C“ = 


0 


0 


-1 


0 


1 


0 


(29) 






0 


0 


0 


-1 


-1 


1 








-1 


0 


0 


0 


0 


-1 





where, in general, C“ has b columns and tit rows. It is apparent that for each column, 
which is associated with a node, 

E'’°li9* = 0 (30) 

k 

or more formally 



C“Q = 0 (31) 

consistent with flow rate conservation at all nodes (see eq (26)). In order to reduce 
the order of the problem by 1 it is useful to define one of the nodes as the ground or 
datum reference. In principle it should not matter which node is taken as this datum 
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node, however, in the coded implementation it is assumed that the circuit nodes 
are numbered such that the datum node is the node bearing the highest numerical 
designator. For example, as seen on Figs 4a, b, the node-5 is the datum node. Because 
the assignment of nodes is independent of the assignment of branches and is arbitrary, 
there is no sacrifice in the generality of the method. Following this point, the node- 
branch incidence matrix, C, can be defined according to the rules dictated by eq (27) 
with the index i satisfies the condition i = 1,2, . . . ,n where 

n = nj — 1 (32) 



Thus, for this example, the node-branch incidence matrix is 

■ 1 1 0 0 0 0 ' 

0-1 1 1 00 

^“0 0-1 0 10 

0 0 0 -1 -1 1 . 

In general C has b columns and n rows. It follows from eq (30) that 



(33) 



CQ = 0 



(34) 



which forces a condition of flow rate conservation or continuity at all nodes. 

The node pressures, defined with respect to the reference datum node, can be 
represented by the column vector 

' Pi 

P= 

. Pn 

which, for the example under consideration, the node pressure vector P is a 4 x 1 
column vector since n = 4. The node branch incidence matrix can be used to predict 
the head losses from the rule 

H = C^P (36) 



(35) 
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where the branch head loss vector, represented as: 



H = 



■ hi ■ 

h2 

. hs . 



( 37 ) 



can be calculated using eqs (11) and (12). For the example being examined, H is a 
6x1 column vector, corresponding to the number of branches. 

In Section 4 the pressure drop vector was introduced without any assumptions 
regarding the physical modehng introduced in Section 2. Under the practical as- 
sumptions made leading to eq (6), it is correct to let 



AP — 



in eq (25), which leads to 



Q = Q,-^YH 



Premultiplication of this equation by C leads to: 



(38) 



(39) 



0 = CQ, + CYH (40) 

where the left hand side of eq (40) is predicted from eq (34). After substitution of eq 
(36) into eq (40) and solving for P, it follows that 

p = y;^q (41) 



where 



Q = -CQ, 



and 



Y„ = CYC^ 



(42a) 



(42b) 
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Given the matrix Y„ and the source vector Q,, any number of matrix inverting 
schemes [Hamming (1962)] could be employed to evaluate eq (40). As suggested in the 
Elhs (1993) thesis the very efficient Cholesky reduction algorithm is applicable here 
because the matrix Y„ is, not only symmetric, but positive definite. The Cholesky 
reduction into upper and lower triangular matrices permits a rapid matrix inversion 
using Gaussian elimination as described in Appendix A. 

Once the node pressures given by eq (40) are obtained the branch head losses and 
flow rates are easily calculated from eq (36) and eq (39), respectively. Note that, once 
the flow rate vector is obtained, the flow rate velocities can be computed according 
to the rule 






V = 




(43) 



vb J 



in agreement with eq (1). The corresponding velocity magnitudes, (|i>i|, |u 2 |, • • • , I^^bI) 
can then be applied with eq (4) to calculate the Reynolds numbers, the friction factors, 
and the turbulence factors using eq (14). The pipe resistances are calculated from eq 
(15) and the corresponding admittances, • • ^Vb, follow from 

1/ri 0 0 ••• 0 



Y = 



0 

0 



l/rj 



0 

0 

0 



(44) 



.0 0 •••0 1/rtj 

If the current turbulence factors differ significantly from the previous turbulence fac- 
tors according to the rule 



MAXERR > max |V>jt(previous) — V’fc(current)| ; {branches k = 1,6} (45) 



the process is repeated. The Fortran input parameter MAXERR defines the maximum 
deviation in the turbulence factors between iterations. It is assumed initially, in the 
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first step of the iterative scheme, that the flow is in the laminar regime, that is, all 
turbulence factors = 1. 

6. Input/Output: Parameters, Conventions, 
Procedures 

It follows for the definition of specific weight that 

1 = P9 (46) 

where g, the gravitational constant, is 32 ft/sec^. Using eq (46) the computer program 
computes the density which is needed for calculation of the Reynolds numbers eq (8). 

Figure 5 defines the convention for distinguishing the entry (F) and exit (T) nodes 
of a particular branch. The cirrow defines the direction of positive flow. 



F T 

• ^ 

Figure 5: Convention defining nodes. 

The required inputs to the FORTRAN program code flow. for provided in Ap- 
pendix A are displayed in Table 1 . 

The calculations axe done in formal English imits of (ft. Lb, sec). However, the 
inputs and outputs are, when appropriate, given in conventional units. Conversion 
factors for these cases are shown in Table 2. 
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Table 1: Input Variables 



Symbol 


Section 


Units 




Reo 


§3 





typical 3500 


MAXERR 


§5 


— 


typical 0.05 


HLCF 


§2 


— 


typical 60.0 


Tlx 


§5 


— 


depends on network 


b 


§5 


— 


depends on network 




§2 


10“® Ib-sec/ft^ 


typical H 2 O, 3.75 


w 


§6 


Ib/fC 


typical H 2 O, 64.0 


node specification —F, —T 


§6 


— 


depends on pipe 


L 


§2 


ft 


depends on pipe 


d 


§2 


inches 


depends on pipe 


e 


§2 


mil 


depends on pipe 


source specification Q, 


§5 


gal/min 


depends on pump(s) 



Table 2: Conversion Factors 



Description 


Symbol 


Formal 

Units 


Conventional 

Units 


Conversion 

Factor 


Flow rate 


Q 


ft^/sec 


gal/min 


450 


Pressure 


P 


Ib/ft^ 


psi 


1/144 


Roughness 


e 


ft 


mil 


12,000 



In order to facilitate the process of batch input, the code has been designed to 
read the input file flowinp.dat rather than require the user to laboriously submit 
the data in the interactive mode. In this way a maximum amount of flexibihty can 
be incorporated into the program without creating a tedious data entry procedure for 
the user of the code. Thus, the “user friendly” goal has been kept in mind. 
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The outputs from the code flow. for, which are listed by branch number (BR), 
are displayed in Table 3. 

Where appropriate both text symbol ajid fortran symbol have been specified in 
Table 4. For example, PF and PT refer to node pressures at entry and exit, respec- 
tively. In some practical cases it may be desirable to employ the equivalent of an 
“ideal” constant flow rate pump, that is, qk = q»k in eq (24). This is readily handled 
by forcing yk 0 or equivalently rjt — ^ oo. Because of the head loss coupling factor 
rule-of-thumb introduced at the end of Section 3, the most convenient way to force 
the pump source to behave ideally is to let dk —* 0. In practice, because of numerical 
instabilities, if yk = 0, it is recommended that the user set 

dfc = 10“^ X min{dj} t = l,2, ...,6 buti/I:. (47) 

A hardcopy of the input/output for the piping networks tested is provided in 
Appendix B. A more detailed description of these cases is covered in the next section. 



Table 3: Output Variable 



Symbol 


Section 


Units 


BR 


§5 


_ 


Q 


§5 


gal/min 


PF 


§6 


10“^ psi (mpsi) 


PT 


§6 


10“^ psi (mpsi) 


HEAD = H 


§5 


10“^ psi (mpsi) 


ERR 


see RHS eq (44) 


— 


Re 


§2 


— 


TBF = V- 


§2 


— 


Y 


§5 


(gaJ/min)/mpsi 
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7. Examples 



General Observations 

Several general observations can be easily confirmed from the data (see Appendix B 
for all the case examples discussed in this section). Specifically, 

PF - PT = HEAD (pf ~Pt = hk) 

and 

Q = Qj + Y HEAD {qk = Qsk + Vkhk) 

In all the examples the sink node of the pump was arbitrarily taken as the datum 
node for convenience of interpretation. The node pressures under columns PF and 
PT are as expected, strictly positive. In all three examples some or all of the pipeline 
segments were operating in the nonlaminar regime. This is evident by checking the 
TBF columns of each example in Appendix B. The associated Reynolds number is 
also provided. The branch errors were calculated according to eq (45). It is easily 
checked that the MAXERR parameter, fine 2 of the input set, put a ceiling on these 
branch error values. 

CASE A — A four branch pipeline network with an ideal source (see Fig 6a and 6b). 
This case illustrates the generation of an ideal flow source. Note, a5 discussed, this 
condition can be obtained by making the pipe diameter of the source br 2 inch relatively 
smedl. The associated Reynolds number and turbulence factors will be artificiaJly high 
and may not, as represented in the fine 1 output for Case A, be within the defined 
format. 



19 



I 2 




CASE B — A six branch pipeline network with nonideal source (sefe Fig 4a and 4b). 
Note, as seen here, the effect of the nonzero.admittance of node 1 is to impede the flow 
of the source into the rest of the pipeline circuit. Also, the source is required to be 
negative because the pump is physically driving the flow opposite to the convention 
defined by the ^Lrrow of branch #1. 
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CASE C — A 12 branch piping network for shoebox coohng rack (see Fig 7). Note 
that several of the branches have negative Q values. The physical interpretation is 
that these cases have a fluid flow opposite to the convention defined on Figure 7. 




Figure 7: 12 branch shoebox coohng rack. 



8. Fiiture Enhancements in the Modeling 

Although the work described herein extends the Ellis (1993) thesis by including the 
nonlaminar flow regime, it still should be considered as preliminary. For example: 

• The network formahsm discussed herein should be compared" with the bench- 
mark Hardy Cross method [Lindeburgjl987)] as well as other more recently 
proposed schemes [CarnaJian and Christensen (1972), Bending and Hutchinson 
(1973), Gay and Preece (1975, 1977)] in order to assess the relative computa- 
tional efficiencies. 

• The next major modehng development in the algorithm would be the addition 
of the capability of evaluating the node to node heat transfers. 
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• Using the theory for compressible fluid flow [Vennard (1961)] required, for ex- 
ample, when the cooling fluid is a gas, the developed scheme could be extended 
to allow for compressible fluids. Due to the sensitivity of gais properties to tem- 
perature this development would have to include effects due to heat transfer. 

• An alternate but similar mathematical formahsm could be derived for pumps 
best characterized ^ls constant pressure sources. It is certainly possible to 
“Theveninize” the source in the context of the present mathematical formal- 
ism. However, it should be noted that because of the nonhnearity inherent in 
the model, this Theveninization would need to be repeated iteratively, decreas- 
ing the efficiency of the process. 

• As previously noted in Section 3, a significant improvement in the rate of con- 
vergence of the method was obtained by reducing the discontinuity in the fric- 
tion factor derivatives at the transition edges. It should be possible to apply 
mathematical adjustments on eq (19), for example, using a technique known as 
Hermite or “osculating” interpolation [Hamming (1962)], in order to completely 
eliminate the slope discontinuity. The potential gain in the rate of convergence 
could be significant. 

• Various mundane, but no doubt practical embellishments could also be consid- 
ered. For example, an option of metric system units could be added. Also, an 
option to specify pipe bends, vertical elevation, and cross-section types could be 
added. If commercial viability is of concern, the program could be dovetailed 
with CAD software. 

• The required background investigation revealed that despite a long history of 
developments on modeling fluid flow in pipes, only recently have mathematical 
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methods involving the theory of chaos begun to unravel and predict the prop- 
erties of nonlaminar fluid flow. A small-scale physical model which predicts all 
the measured properties on the Stanton diagram is at this date unavailable. 

9. Conclusion 

A mathematical model for predicting laminar and nonlaminar flows in pipehne net- 
works hcis been proposed and successfully tested. This scheme, upon iterative apph- 
cation, has been found to converge to an approximate solution. The computational 
efficiency of this convergence process was found to be highly dependent on the fric- 
tion factor modeling in the transition region between laminar and turbulent behavior. 
In particular, a significant improvement in the convergence rate was obtained by re- 
placing the previously proposed linear model with a nonlinear model having a less 
dramatic first derivative discontinuity. A linear gauge for the deviation from laminar 
behavior, referred to as the turbulence factor, was introduced in this report in order 
to facilitate the computational task. 
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10: Nomenclature 



Roman Letter Symbols 

A pipe cross-sectional area 

d pipe diameter 

e pipe roughness 

/ friction factor 
g gravitational acceleration 
h head loss due to friction 
P pressure 
Q flow rate 

R pipe resistance 

Re Reynolds number 

V velocity 

Ws time-rate-change of work output (pump) 

Y pipe admittance 
z elevation 

Greek Letter Symbols 

7 specific weight 

e pipe roughness ratio 

/X viscosity 
p density 



Units 

ft^ 

ft 

ft 

ft / sec^ 
lb/ft2 
lb/ft2 
ft^/ sec 
Lb- sec/ft® 

ft /sec 
Ib-ft/sec 
ft^/(lb-sec) 
ft 



lb/ft3 

Ib-sec/ft^ 

slug/ft^ 
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Appendix A: Cholesky Reduction 



This appendix provides a more complete description on the Cholesky reduction algo- 
rithm [Kraus (1992)]. Positive definite, symmetric matrices can be factored according 
to the rule 

M = LL^ (Al) 

where L is an upper triangular matrix. Because of the applicabiUty of this theorem 
[Ellis (1993), Kraus (1987)] to the network matrix Y„ eq (41b) it follows that 

Yn = LL^ (A2) 



and for eq (40) 



LL^P = Q 



(A3) 



The advantage of the factorization eq (Al) is now apparent since eq (A3) can be 
broken into two parts: First, 

LJ = Q (A4a) 



and second, 

L^P = J (A4b) 

which are easily solved in succession via Gaussian elimination to determine first the 
column vector J as an intermediate step then the derived node pressures P. 
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Appendix B: Input/Output Hardcopy 



**«**«*««**«**«* INPUT FORMAT o************************ 

REYNOLDS’ NUMBER FOR RIGHT EDGE OF TRANSITION REGION 

ERROR TOLERANCE IN CALCULATION 

THE HEADLOSS COUPLING FACTOR (TYPICAL NUMBER 60 ) 

THE NUMBER OF NODES 
THE NUMBER OF BRANCHES 
THE VISCOSITY X 10 -5 LB-SEC/FT'3 
THE DENSITY LB/FT ‘3 

STARTING NODE, ENDING NODE, LENGTH(FT) ,DIAMETER(IN) , ROUGHNESS (MILS) 

NOTE FOR IDEAL CURRENT SOURCE SET DIAM=.001X SMALLEST 
THE NUMBER OF SOURCES 
SOURCE STRENGTH GAL/MIN 
SOURCE BRANCH 

*nnnnLif^Lif>if*************************************************** 

INPUT (CASE A 4 BRANCHES IDEAL SOURCE) 

3500.0 
0.005 

60.0 
4 

4 

3.75 

64. 

4.1, 1.0, .001,1.0 

1 . 2 , 1 . 0 , 1 . 0 , 2.0 

2.3, 1.0, 1.0, 0.5 

3.4, 1.0, 1.0,0. 6 
1 

6.0 

1 



OUTPUT CASE A 



6R 


Q 


PF 


PT 


HEAD 


ERR 


Re 


TBF 


Y 


1 


6.00 


.0 


282.1 


-282.1 


.0039 


0i*t*t** 


4c 4c * 4c 4c 4c 


.0000 


2 


6.00 


282.1 


183.8 


98.2 


.0032 


10806.2 


5.6 


.0609 


3 


6.00 


183.8 


92.1 


91.7 


.0032 


10806.2 


5.2 


.0652 


4 


6.00 


92.1 


.0 


92.1 


.0032 


10806.2 


5.3 


.0649 



***************^ 1 :^ 1 ****************************************************** 
INPUT (CASE B 6 BRANCH PIPELINE NETWORK) 

3500.0 
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0.01 

60.0 

5 

6 

3.75 

64. 

1,5, 1.0, 1.0, 1.0 

1 , 2 , 1 . 0 , 1 . 0 , 2.0 

2.3, 1.0, 1.0, 0.5 

2.4, 1.0, 1.0, 0.6 

3.4, 1.0, 1.0, 0.6 

4.5, 1.0, 1.0, 5.0 
1 

- 6.0 

1 

OUTPUT CASE B 



BR 


Q 


PF 


PT 


HEAD 


ERR 


Re 


TBF 


Y 


1 


-1.75 


25.1 


.0 


25.1 


.0073 


3149.7 


2.0 


.1679 


2 


1.75 


25.1 


14.6 


10.6 


.0075 


3149.7 


2.1 


.1645 


3 


.58 


14.6 


12.9 


1.7 


.0000 


1049.9 


1.0 


.3409 


4 


1.17 


14.6 


11.2 


3.4 


.0000 


2099.8 


1.0 


.3409 


5 


.58 


12.9 


11.2 


1.7 


.0000 


1049.9 


1.0 


.3409 


6 


1.75 


11.2 


.0 


11.2 


.0078 


3149.7 


2.2 


.1554 



CASE C ( 12 BRANCH SHOEBOX COOLING RACK) 

3500.0 

0.01 

60.0 
8 

12 

3.75 

64. 

8 , 2 , 1 . 0 , 1 . 0 , 1.0 

2.3, 1.0, 1.0, 2.0 

3.4, 1.0, 1.0, 0.5 

4.8, 1.0, 1.0, 0.6 

5.8, 1.0, 1.0, 0.6 

2.6, 1.0, 1.0, 5.0 

7.3, 1.0, 1.0, 1.0 

1.4, 1.0, 1.0, 2.0 

6.5, 1.0, 1.0, 0.5 

6.7, 1.0, 1.0, 0.6 
7,1, 1.0, 1.0, 0.6 
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5,1, 1.0, 1.0, 5.0 
1 

6.0 

1 

OUTPUT CASE C 



BR 


Q 


PF 


PT 


HEAD 


ERR 


Re 


TBF 


Y 


1 


3.48 


.0 


24.6 


-24.6 


.0341 


6273.0 


3.4 


.0990 


2 


1.76 


24.6 


14.5 


10.1 


.0393 


3168.5 


2.0 


.1681 


3 


1.33 


14.5 


9.6 


4.9 


.0084 


2394.5 


1.3 


.2685 


4 


1.74 


9.6 


.0 


9.6 


.0371 


3138.2 


1.9 


.1750 


5 


1.74 


9.6 


.0 


9.6 


.0366 


3134.8 


1.9 


.1752 


6 


1.72 


24.6 


14.4 


10.2 


.0386 


3104.5 


2.1 


.1635 


7 


-.43 


13.2 


14.5 


-1.3 


.0000 


773.9 


1.0 


.3409 


8 


.41 


10.8 


9.6 


1.2 


.0000 


743.6 


1.0 


.3409 


9 


1.32 


14.4 


9.6 


4.8 


.0072 


2381.3 


1.3 


.2710 


10 


.40 


14.4 


13.2 


1.2 


.0000 


723.2 


1.0 


.3409 


11 


.83 


13.2 


10.8 


2.4 


.0000 


1497 . 1 


1.0 


.3409 


12 


-.42 


9.6 


10.8 


-1.2 


.0000 


753.5 


1.0 


.3409 
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Appendix C: Computer Code 



$debug 

$LIST 

c program for CALCULATING PRESSURES AND FLOW RATES IN PIPILINE NETWORKS 



C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 



********** ************************** ***^****** ««****«* ************* 

PROGRAM flow. FOR WRITTEN BY PROF. RON J PIEPER 

NPS 408 6562101 
AUG. 30, 1994 

SEE FILE flowINP.DAT FOR INPUT 

SEE FILE flowOUT.DAT FOR OUTPUT 

SEE FILE flowDIA.DAT FOR DIAGNOSTICS ON INPUT 

DESIGNED TO HANDLE UP TO 40 BRANCHES 

CAN BE ADJUSTED BY INCREASING DIMENSION 
OF THE ARRAYS 

******************************************************************** 



c 


mu 


fluid viscosity (lb-sec/ft*2) 


c 


rho 


fluid density (lb/ft*3) 


c 


g 


acceleration due to gravity (ft/s 


c 


nu 


specific gravity rhoxg 


c 


pi 


pi 


c 


ren# 


reynold’s t (rho*vel*d/mu) 


c 


f 


friction factor 


c 


vel 


velocity of fluid 


c 


pres 


pressure of the fluid 



c 

c 

c 

c 

c 



*****************^i^e veiriables ************************************** 



ell 

d 

e 

rat 



pipe lengtb(s) 
pipe diameter(s) 
roughness 
roughness ratio 



C HCLF 

Q ********** 

C TURB 
C MAXERR 
c 
C 
c 



HEADLOSS COUPLING FACTOR 
F FACTOR **♦*♦*♦♦**♦♦♦♦♦♦♦*********♦**♦*♦ 

TURBULENCE FACTOR = F*2100/64 (=1.0 LAMINAR REG) 
MAXERR IN TURBULENCE FACTOR 

FOR REFERENCE ALSO SEE ITERATIVE COOLBROOKE-WHITE 
AND ALSO SEE DIRECT CHURCHILL eQ ( IN PLOT PROGRAM) 
applies for rent > REYEDGE 



C REDGE THE LEFT EDGE OF TRANSITION REGION 

C REYNOLDS# WHERE COOLBROOK BEGINS 



c ************************************************************** 

c *i»***i^ network constants ************************************* 
c N the nufflver of junctions in the system 
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c N1 
c B 
C S 



the # junction - datum node (n-1) 
the i branches in the system 
the # of sources 



C :^******J¥***************************************************^ 



C **********************************t***********0*t0*t ********* 



integer N,N1,B,QB,S,DN 
integer NT(40), NF(40),X, itest 

real mu, rho, g, pi, C(40,40), d(40), ell(40), r(40) ,QS1 ,QS(40, 1) 
REAL RLAM(40),TURB(40),Ct(40,40),Cy(40,40),YLAM(40), Y(40,40) 

REAL IS(40,1),YN(40,40),P(40,1), YH(40,1) ,V(40,1) , REY(40,1) 

REAL AREA(40) ,ISI(40,1) ,YNI(40,40) ,Q(40,1) ,ERR(40) ,ep(40) ,rat(40) 
REAL TURBF , MAXERR , QG (40 ) , FCROS (40 ) 

REAL REDGE, EPS,H(40, 1) ,HCLF 
CHARACTER*20 FORM 
CHARACTER*3 BR 
CHARACTER*? NAME(ll) 

BR=’BR ’ 

NAME(7)=’ TBF ’ 

NAME(4)=’ HEAD' 

NAME(2)*' PF ’ 

NAME(3)*' PT ’ 

NAME (6)=’ Re ’ 

NAME(1)=> Q’ 

NAME(8)=’ Y ' 

NAME(5)=’ ERR’ 

4 FORMAT (A12) 

OPEN (UNIT*?, FILE* ’flowDIA.DAT’ , STATUS* 'NEW' , ACCESS* ’ SEQUENTIAL ’ ) 
OPEN (UNIT*8 , FILE* ’ f lowINP . DAT ’ , STATUS* ’ OLD ’ , ACCESS* ’ SEQUENTIAL ’ ) 
OPEN (UNIT*9, FILE* 'flowOUT.DAT', STATUS* ’ NEW ’, ACCESS* 'SEQUENTIAL') 

0t************************************************************** 

g*32.1?4 
pi*3. 1415926 
READ (8, 8) 

8 FORMAT 

C write(*,13) 

write(?, 13) 

13 format ('INPUT THE REYNOLDS i FOR RIGHT EDGE TRANSITION ’) 



c 

C 

C 

c 



c 

c 



ASSUMPTIONS 

1. COOLBROOKE-WHITE FORMULAE FOR REY#> redge 

2. WORK IN ENGLISH UNITS 

3. FLUID IS INCOMPRESSIBLE 

4. ALL PIPES ARE CIRCULAR IN SHAPE 

5. THE CORRECT SOLUTION IS OBTAINED BY ITERATION 
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read(8,30) redge 
WRITEC7.30) REDGE 
C write(*,23) 

¥rite(7,23) 

23 format (' enter the maximum error IN TURBULENCE FACTOR desired’ ) 

read(8,30) MAXERR 
WRITE (7, 30) MAXERR 
C WRITE(*,21) 

21 FORMAT (’ENTER THE AVERAGE HEADLOSS COUPLING FACTOR, TYPICAL #60’) 

READ (8, 30) HLCF 
write(7,21) 
write (7, 30) HLCF 
C write(*,10) 

write(7,10) 

10 format (’ enter the number of nodes in the rack system N>100 ’ ) 
read(8,15) N 
N1=N-1 

15 format (15) 

urite(7,15) N 
C write(*,15) N 

C write(*,20) 

write(7,20) 

20 format (’ enter the number of branches in the rack system B>N ’ ) 
read(8,15) B 
C write(*,15) B 

write(7,15) B 

********t************************************************************** 

C write(*,25) 

write(7,25) 

25 format (’ input the viscosity X 10“-5 lb-SEC/ft*2 ’ ) 

read(8,30) mu 
WRITE (7, 30) mu 
C WRITE(*,30) mu 

mu*mu* .00001 
30 format (flO.4) 

C write(*,35) 

write(7,35) 

35 format ( ’ input the specific weight lb/ft‘3’ ) 
read(8,30) sw 
rho=sw/g 

C write(*,30) sw 

write(7,30) sw 
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c 



Node-branch matrix 



c 

c matrix initialization 
c 

do 40 i=l,Nl 
do 40 j=l,B 
40 C(i,j)=0 

do 42 i=l,B 
NF(i)=0 
42 NT(i)=0 

C 
C 

C receive dat from the keybord to develop a,d and ell matrices 
c for circuleu: pas ages 

c start loop on branches to input node to node information 

do 50 i=l,B 
Hrite(7,56) i 

55 format (* at branch number 2x, i5) 

C write(*,60) 

WRLTE(7,60) 

60 formatC'ENTER nf (i) ,nt(i) ,length(ft) .diameter(in) ,ROUGHNESS-mil’) 

read(8,*) NF(i) ,NT(i) ,ell(i) ,d(i) ,ep(i) 
ep(i)=ep(i)*.001 
rat(i)=ep(i)/d(i) 
c convert from mils to inches 

WRITE(7,*) NF(i) ,NT(i) ,ell(i) ,d(i) ,ep(i) ,rat(i) 

C WRITEC*,*) NF(i),NT(i),ell(i),d(i),ep(i),rat(i) 

C 70 formate i5,i5,f8.3,f8.3,f9.5) 

c ************* calculate the coolbrook crossing value at R= 3500 
c colebrook-Bhite formula page 198 of fox, let f goto 4f 

sqf= (64. 0/2100.0) **.5 
3100 save*sqf 

sqf=l .0/(-2.0*logl0(2.51/(REDGE*save) + rat(i)/3.7)) 

IF (SAVE .EQ. SQF) GOTO 4000 
error® abs (save-sqf) /sqf 
if (error .gt. .001) goto 3100 
4000 f cros(I)=sqf**2 

C WRITE(*,*) >I,FCR0S(I),RAT(I)’, I. FCROS(I), RAT(I) 

50 continue 

C SECTION ON SOURCE INPUT 
c INPUT SOURCES, S OF THEM 

C QS MATRIX OF FLOW SOURCES 

C QB SOURCE BRANCH 
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C QSl SOURCE STRENGTH 

C WRITE(*,90) 

WRITE (7, 90) 

90 FORMAT ( ' INPUT THE NUMBER OF SOURCES ') 

READ(8,92) S 
write (7, 92) S 
C write(>»,92) S 

92 FORMATdS) 

DO 98 1*1, S 
C WRITE(*,95) I 
WRITE(7,95) I 

95 FORMAT( ’INPUT THE SOURCE STRENGTH FOR THE’, 15, ’TH SOURCE’) 
READ(8,30) QSl 
write(7,30) QSl 
C write(*,30) QSl 

QSl=QSl/450 

C CONVERT FROM GAL/MIN TO FT*3/SEC CONVERSION FACTOR 1/450 
C X GAL/MIN* (450)“-! FT*3/SEC 
. WRITE(7,93) 

C WRITE(*,93) 

93 FORMAT (’ INPUT THE BRANCH NUMBER OF THE SOURCE ’) 

READ(8,15) QB 
write(7,15) QB 
C write(’»,15) QB 

98 CONTINUE 

DO 100 K=1,B 
IF (K .EQ. QB) THEN 
QS(K,1)= QSl 
ELSE 

QS(K,1)=0 
ENDIF 
100 CONTINUE 

c Staart setting up ’ C ’ matrix 
do 120 i*l,B 

if (NF(i) .LT. N) then 
C(NF(i),i)*l 

endif 

if( NT(i) .LT. N) then 
C(NT(i),i)=-l 

endif 

c convert diameter in inches to feet 
d(i)*d(i)/12.0 



33 



area(i)=pi*(d(i)**2)/4 

length=ell(i) 

c 60 diameters accounts for in/out head loss 
addl=HLCF*d(i) 
ell(i)=length+ addl 

C CALCULATE THE LAMINAR RESISTANCES FOR EACH BRANCH 
RLAH(i)= 32.0* mu* ell(i)/(AREA(I) * (d(i))**2) 
YLAH(I)=1/RLAM(I) 

120 CONTINUE 

C WRITE "C" MATRIX TO SCREEN AND RECORD FOR CHECK 
WRITE (7, 123) 

C WRITE (*,123) 

123 FORMAT ( ' THE C MATRIX BELOW. TRANSPOSED ’ ) 

DO 126 1=1, B 

C WRITE(*,125) (C(J,I),J=1,N1) 

WRITE(7,125) (C(J,I),J=1,N1) 

125 F0RMAT(8(1X,F8.4)) 

126 CONTINUE 

C WRITE THE LAMINAR RESISTANCE AND ADMITTANCE VECTORS 

WRITE(7,*) ’THE LAMINAR RESISTANCE AND ADMITTANCE VECTORS’ 

C WRITE (*,*) ’THE LAMINAR RESISTANCE AND ADMITTANCE VECTORS’ 

DO 130 1=1, B 

WRITE(7,*) RLAM(I),YLAM(I) 

C WRITE(*,*) RLAM(I) ,YLAM(I) 

TURB(I)=1.0 

130 CONTINUE 

C MATRIX MULTIPLICATION OF MATRICES 

C 

C 

DO 140 J=1,B 
DO 140 K=1,B 
Y(J,K)=0 

IF (J .EQ. K) THEN 

Y(J,J) = YLAM(J) 

ENDIF 

140 CONTINUE 

C MATRIX MANIPULATION SECTION 
c N1 rows of C 
c b columns of C 
c CT transpose of C 

CALL TRANSP(C,CT,N1,B) 

C write(*,*) ’ write transpose of matrix C ’ 
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write(7,*) ’ write transpose of matrix C ’ 
call Warray(CT.B.Nl) 

C srite(*,*) ’ write matrix Y’ 
write (7,*) ’ write matrix Y’ 
call Warray(Y.B.B) 



C PROPOSED ENTRY POINT FOR ITERATIVE SCHEME 
250 call matmul(Cy,C,Y,Nl,B,B) 

C writeC*,*) ‘ write PRODUCT C*Y NIXB ’ 
write(7,*) ’ write PRODUCT C*Y NIXB ’ 
call Warray(Cy,Nl ,B) 

CALL HATMUL(YN,Cy,Ct,Nl,B,Nl) 
c node flow source vector CQs 

call matmuK Is ,C,QS,N1 ,B, 1) 
c commented out ao-tifact of old scheme 
DO 200 1=1, N1 
IS(I,1)=-IS(I,1) 

ISI(I,1)=IS(I,1) 

200 CONTINUE 

DO 295 1=1, N1 
DO 295 J=1,N1 

YNI(I,J)=YN(I,J) 

295 CONTINUE 

C SECTION ALSO WRITES THE ARRAY PRIOR TO BEING cHOLESKY PROCESSED 
WRITE (7,*) ’YNI MATRIX PRIOR TO CHOLESKY * 

C WRITE (*,♦) ’YNI MATRIX PRIOR TO CHOLESKY ’ 

call Warray(YNI,Nl,Nl) 

WRITE (7,*) ’ISI MATRIX PRIOR TO CHOLESKY » 

C WRITE(*,*) ’ISI MATRIX PRIOR TO CHOLESKY ’ 

call Warray(ISI,Nl,l) 
call CH0LESKY(YNI,ISI,N1) 

WRITE (7,*) ’YNI MATRIX FOLLOWING CHOLESKY ’ 

C WRITE(*,*) ’YNI MATRIX FOLLOWING CHOLESKY ’ 

call Warray(YNI,Nl,Nl) 

C P NODE MATRIX P=YN“-1*IS 
C H BRANCH PRESSURE CT*P 
C Q BRANCH FLOW RATE 
C YH MATRIX PRODUCT Y*H 
DO 296 1=1, N1 
P(I,1)=ISI(I,1) 

296 CONTINUE 
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CALL KATMUL(H,CT,P,B,N1,1) 

CALL MATMUL(YH,Y,H,B,B,1) 

DO 300 1=1, B 
Q(I,1)=YH(I,1)+QS(I,1) 

C DIVIDE PRESSURES IN LB/FT* 2 BY 144 TO CONVERT TO PSI 
300 CONTINUE 

c REYNOLDS NUMBER REGIME DETERMINATION 
DO 5000 1=1, B 

V(I,1)=ABS( Q(I,1)/AREA(D) 

REY(I , l)=rho*V(I , l)*d(I)/mu 
REN=REY(I,1) 

C WRITE(*,*) ’ REYNOLDS REY(I,1) 

WRITE(7,*) ' REYNOLDS #', REY(I,1) 
if (REY(I,1) .LE. 2100.0) THEN 
Y(I,I)= YLAM(I) 

TURB(I)»1.0 
TURBAV=1.0 
ERR(I)=0 
F=64. 0/2100.0 
else 

if (REY(I,1) .IE. REDGE ) THEN 
X1=REDGE-2100.0 
Yl=Fcros(i) - 64.0/2100.0 

F= 64.0/2100.0 ♦ Y1*SIN(PI* (REN-2100)/(2*X1)) 
else 

c colebrook-white formula 

sqf=(64. 0/2100. 0)**. 5 
1310 save=sqf 

8qf=1.0/(-2.0*logl0(2.51/(Ren*save) + rat(i)/3.7)) 

C IF (SAVE .EQ. SQF) GOTO 1400 

error=abs (save-sqf ) /sqf 
if (error .gt. .001) goto 1310 
1400 f=sqf**2 

endif 

TURBF=F*Ren/64.0 
TURBAV= (TURB (I ) +TURBF) /2 . 0 
Y ( I , I ) = 1 / (RLAM ( I ) *TURBAV) 

ERR ( I ) *ABS (TURB AV-TURB ( I ) ) /turb ( i ) 

TURB(I)=TURBAV 

ENDIF 

C WRITE(*,*) I, 'FRICTION FACTOR’, F,’THE TURB FACTOR’ , TURB AV 
WRITE(7,*) I, ’FRICTION FACTOR’, F,’THE TURB FACTOR’ .TURBAV 
5000 CONTINUE 
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EPS=0 

DO 5100 J=1,B 

C WRITE(*,*) ’ THE ERROR FOR BRANCH’ ,J,’ IS’ ,ERR(J) 

WRITE(7,*) ’ THE ERROR FOR BRANCH’, J.’ IS’.ERR(J) 

IF (ERR(J) .GT. EPS ) EPS=ERR(J) 

5100 CONTINUE 

IF (EPS .GT. MAXERR) GOTO 250 

C C CCCCCCC NOTE ERR DEFINED IN TERMS OF TURBF IS ALREADY NORMALIZED 
WRITE(9,5900) BR.NAHE(l) ,NAME(2) ,NAME(3) ,NAME(4) , 

♦ NAME(5),NAME(6),NAHE(7),NAME(8) 

5900 F0RMAT(A3,8(1X,A7)) 

C 

C READY TO BEGIN OUTPUT OF SOLUTION 

C WRITEC*,*) ’BRANCH, FLOW(Q).PF(MPSI),PT(MPSI), head, ERR, 

C ♦ REYi, TURB FACTOR, FLOW(Q-QS) , ’ 

WRITEC7,*) ’ BRANCH, FLOH(Q), presURE F, pres T, head, ERR, 

+ REY#, TURB FACTOR, FLOH(Q-Qs)’ 

DO 7000 J=1,B 

c Sources calculated in terms of effective pressure diff ( thevenin) 
C NOTE THE FLOW VALUES ARE IN GAL/MIN 
C CONVERT FLOW VALUES TO GAL/MIN 

C DIVIDE PRESSURES IN LB/FT*2 BY 144 TO CONVERT TO PSI 
C ADMITTANCE Y CONVERTED TO (GAL/MI N)/mPSI 
c * 1000 to convert to mpsi 

PF=P(NF(J),1)/144.0*10**3 

PT=P(NT(J),1)/144.0'»10**3 

H(j,l)=H(j,l)/144.0*10**3 

Y(J,J)=Y(J,J)*450.0*144.0/10**3 

QS(J,1)=450.0*QS(J,1) 

Q(J,1)=Q(J,1)*450 
C QG(J)=Q(J,l)-Qs(J,l) 

C HRITE(*,6000)J,Q(J,l),PF,PT,H(j,l),ERR(j), 

C ♦ REY(J,1),TURB(J),Y(J,J) 

HRITE(7,6000)J,Q(J,l),PF,PT,H(j,l),ERR(j), 

♦ REY(J, 1), TURBO) ,Y(J,J) 

HRITE(9,6000)J,Q(J,l),PF,PT,H(j,l),ERR(j), 

♦ REY(J,1), TURBO) ,Y(J.J) 

6000 F0RMAT(I3,1X,F7.2,3(1X,F7.1),1X,F7.4,2(1X,F7.1),1X,F7.4) 

7000 CONTINUE 

HRITE(*,*) ’ ***************************************** 

VmiTEC*,*) ’ ’ 

VmiTE(*,*) ’ SEE FILE f 1 ohOUT.DAT FOR OUTPUT ’ 

WRITEC*,*) ’ SEE FILE flouDIA.DAT FOR CHECKING INPUTS ’ 
VmiTEC*,*) ’ ’ 
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WRITE(*,*) ’ 

STOP 

END 

subroutine cholesky(a,b,n) 

C SOLVES THE PROBLEM AX=B 
C X UNKNOWN N TUPLE 

C A KNOWN RANK 2 MATRIX OF ORDER N 
C B KNOWN N TUPLE ( UPON ENTRY ) 

C B=X ( UPON EXIT) 

real b(40,l) 
real a(40,40) 

write (7,*) ’ in cholesky routine, matrix order n follows ’ 

call Warray(a,n,n) 

call decomp(a.n) 

b(l,l)=b(l,l)/a(l,l) 

do 10 i=2,n 

d=b(i.l) 

il=i-l 

do 5 1=1, il 

5 d=d-a(l,i)*b(l,l) 

10 b(i,l)=d/a(i,i) 

b(n,l)=b(n,l)/a(n,n) 

nl=n-l 

do 30 1=1 ,nl 

k=n-l 

kl=k+l 

do 20 j=kl,n 

20 b(k,l)=b(k,l)-a(k,j)*b(j,l) 

30 b(k,l)=b(k,l)/a(k,k) 

return 
end 

subroutine matmul(c,a,b,n,m,l) 

C SOLVES THE PROBLEM C=A*B 

c C UNKNOWN MATRIX Ixn 

c A KNOWN MATRIX MxN 

c b KNOWN MATRIX IXm 

real a(40,40), b(40,40), c(40,40) 

C n i columns of a and c 
cm i rows of a and columns of b 

cl t rows of b and c 

do 25 i=l,n 

do 25 3=1,1 
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25 

410 

1 

2 

3 

10 

20 

21 

30 

40 

45 



c(i,j)=0 
do 25 k=l,m 

c(i,j)= c(i,j)+ a(i,k)*b(k, j) 

return 

end 

subroutine decomp(a,n) 

real a(40,40) 

i=l 

write (7,*) ’ now in decomp routine ’ 
write (7,*) 'the order of matrix a is’, n 
do 410 ii=l,n 

write(7,35) (a(ii, j j) , j j=l,n) 
continue 

if ( a(l,l) ) 1, 1, 3 
write(7,2) 
write(7,*) i,a(i,i) 

format ( ’ zero or negative radicand ’) 
goto 200 

a(l,l)=sqrt(a(l,D) 

do 10 j=2,n 

a(l, j)=a(l, j)/a(l,l) 

do 40 i=2,n 

il=i-l 

d=a(i,i) 

do 20 1=1, il 

dhold=d 

d=d-a(l,i)*a(l,i) 
if ( d .eq. 0) then 

WRITEC7,*) ’ d is zero ’ 
write(7,*) i,a(l,i), dhold 

endif 

write(7,*) ’ i, matrix(i,i) follow ’ 
write(7,*) i, a(i,i) 
if (a(i,i)) 1,1,21 
a(i,i)*sqrt(d) 

if ( i .eq. n) goto 45 
i2=i+l 

do 40 j»i2,n 
d=a(i, j) 
do 30 1=1, il 
d=d-a(l,i)*a(l, j) 
a(i,j)»d/a(i,i) 
do 50 i=2,n 
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c zero out lower paurt of matrix 

do 50 j=l,il 
50 a(i,j)=0 

write (7,*) ’made it through decomp , matrix follows ’ 

do 400 ii=l,n 

write(7,35) (a(ii, j j) ,j j=l,n) 

35 format(8(2x,f8.4)) 

400 continue 

200 ret\im 

end 

SUBROUTINE TRANSP(A, AT.NR.NC) 

REAL A(40,40), AT(40,40) 

DO 20 J=1,NC 
DO 20 K=1,NR 
20 AT(J,K)=A(K,J) 

RETURN 

END 

subroutine warray(a,n,m) 

c writes array to screen and to diagnostic file 
real a (40, 40) 
c n rows of a 

c m coliimns of a 

do 20 j=l,n 

write(7,30) (a(j ,i) ,i=l,m) 

C write(*,30) (a(j ,i) ,i=l,m) 

20 continue 

30 format (8(lx,f8.4)) 

return 
end 
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